Absence of static stripes in the two-dimensional t—J model determined using an 
accurate and systematic quantum Monte Carlo approach 
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We examine the two-dimensional t—J model by using variational approach combined with well 
established quantum Monte Carlo techniques [S. Sorella et al, Phys. Rev. Lett. 88, 117002 (2002)] 
that are used to improve systematically the accuracy of the variational ansatz. Contrary to recent 
density-matrix renormalization group and projected entangled-pair state calculations [P. Corboz et 
al, Phys. Rev. B 84, 041108(R) (2011)], a uniform phase is found for J/t = 0.4, even when the 
calculation is biased with an ansatz that explicitly contains stripe order. Moreover, in the small 
hole doping regime, i.e., S < 0.1, our results support the coexistence of antiferromagnetism and 
superconductivity. 
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Introduction. The comprehension of the low-energy 
properties of strongly-correlated systems remains one 
of the biggest challenges in modern condensed matter 
physics. Indeed, although a fair good understanding 
has been achieved in some limiting cases (especially for 
large spatial dimensions, thanks to dynamical mean- field 
theory-'-), many important questions remain wide open 
in the two-dimensional case, where the competition be- 
tween charge/spin ordering and superconductivity is very 
strong. Unfortunately, in this case, there are not unbi- 
ased techniques that may be used to obtain accurate re- 
sults for low temperatures and large system sizes. There- 
fore, several approximate methods have been developed 
and applied in the last years, like for example variational 
(VMC)^ and fixed-node (FN) Monte Carlo,"* density- 
matrix renormalization group (DMRG)^ or its develop- 
ments based upon the so-called tensor network states, in- 
cluding multi-scale entanglement renormalization ansatz 
(MERA)^ and projected entangled-pair states (PEPS)'^, 
which has been recently generalized to fermionic systems^ 
and infinite lattices (iPEPS).^ 

Different calculations on the t—J model have shown 
contradicting outcomes^Sr— and the whole phase dia- 
gram of this model is still highly debated. One important 
issue, related to the mechanism of pairing in the cuprate 
materials, is whether some charge instability may take 
place (at q = 0, leading to phase separation, or at finite 
q, leading to the so-called stripes) or instead the homo- 
geneous ground state is stablei^i^ In the latter case, the 
residual attraction among quasi-particles may lead to a 
superconducting state. Previous FN calculations empha- 
sized the existence of a stable superconducting ground 
state^^ while DMRG and iPEPS results suggested a 
stripe order 1^ 

The competition between superconductivity and 
stripes have been studied in several papers and differ- 
ent aspects have been addressed in the recent past^^— 
For example, two of us showed that a relatively small 
anisotropy in the super-exchange (and hopping) param- 
eters may lead to a striped order.^^ In this regard, it is 



crucial to have a controlled method that may give vari- 
ational results, in order to make a direct comparison of 
energies (and other correlation functions) among differ- 
ent methods and reach a final consensus. 

In this Rapid Communication, we adopt the same 
method used in Ref. 20: by applying few Lanczos steps 
to the variational wave function and by filtering out its 
high-energy components (by means of the Green's func- 
tion Monte Carlo with the FN approximation), the ac- 
curacy of the calculations may be highly improved. This 
approach is particularly effective at low doping and is ac- 
tually unbiased at half filling. Moreover, an estimation 
of the exact energy may be given by the variance ex- 
trapolation: besides the energy, also the variance of the 
state can be calculated, and the energy with zero vari- 
ance can be extracted. From our finding, the existence 
of a striped phase for S ~ 1/8 is rather unlikely: even the 
best approximation to the ground state does not show 
any evidence towards charge inhomogeneity. Although 
the present calculations cannot rule out the possibility 
to have small static stripes, our Monte Carlo approach is 
expected to reproduce qualitatively correct ground-state 
properties; in particular, it is reliable for determining 
the spatial dependent hole density: whenever an exter- 
nal modulated potential is added to the t—J Hamilto- 
nian, the FN approximation gives rise to stripes, even 
when the initial state is chosen to be homogeneous. 

Model and methods. The t—J model on the two- 
dimensional square lattice is defined by: 

H^-tY^ 4^cj^ + H.c. + J ^ f S, • - ^n.Tij j , 

(1) 

where (• • • ) indicates nearest-neighbor sites, cj^ (cj,^) 
creates (destroys) an electron with spin a on the site i; 
Si and Ui are the spin and density operators on the site 
i, respectively. The t—J Hamiltonian is defined in the 
subspace without doubly occupied sites. In the following, 
we will take the amplitude for nearest-neighbor hopping 
t = 1, and consider the super-exchange J/t = 0.4. The 



2 



hole doping will be denoted by 5 = 1 — N/L^ where N 
and L are the number of electrons and sites, respectively. 
Periodic boundary conditions are taken in both directions 
and L X L or 45-degree tilted lattices (with L — 2P, I 
being an odd integer, so that the non-interacting ground 
state is non-degenerate at half filling) are considered. 
Our starting variational wave function is defined as 



\-^v) =VNPGJdJs\^MF), 



where Vn is the projector onto the subspace with 
N electrons, Vg is the Gutzwiller projector, which 
enforces no double occupation on each site; J7d = 
exp{l/2J2i^jUijninj) and Js = exp{l/2Y,i^jVijSfS^) 
are density-density and spin-spin Jastrow factors, re- 
spectively. Finally, |$mf} is a mean- field state that 
may contain BCS pairing, antiferromagnetic order, or 
both. In our recent papersji^iiS we have shown that very 
good variational energies can be obtained by orienting 
the magnetic order parameter in the x—y plane, so that 
quantum fluctuations may be included thanks to the Jas- 
trow term J^s- In this case, however, the wave function 
takes the form of a PfaffianJ^ii^ Conversely, whenever 
the antiferromagnetic order is taken along the z direction, 
we deal with a determinant.'^ The variational parameters 
are the Uij 's and Vij 's (for all independent distances in the 
lattice) and few parameters that describe the mean-field 
state |$mf) (i-e., the pairing amplitude Abcs, the an- 
tiferromagnetic parameter Aaf, as well as the chemical 
potential and the next-nearest-neighbor hopping describ- 
ing the variational electron dispersion). Due to the pres- 
ence of strong correlations (i.e., the Gutzwiller projector 
and the Jastrow factors), a variational Monte Carlo ap- 
proach is required to compute the energy and all physical 
observables. 

The accuracy of the wave function Q may be improved 
in different ways. The first one is by applying Lanczos 
steps: 



(3) 



fe=i 



where the a^'s are additional variational parameters. 
Clearly, whenever is not orthogonal to the exact 

ground state, |\l/p) converges to it for large p. However, 
on large sizes, only few steps can be efficiently performed: 
here, we consider the case with p = 1 and p — 2 {p — 
corresponds to the original variational wave function). 
Moreover, an estimation of the true ground-state energy 
may be achieved by the variance extrapolation: for suffi- 
ciently accurate states, we have that E « _Eox+constX(T'^, 
where E = {'H)/L and cr^ ^ ((7^2^ _ {UfVL are the en- 
ergy and variance per site, respectively. Therefore, the 
exact ground-state energy i?cx may be assessed by fitting 
E vs (T^ for p = 0, 1, and 2. 

Another way to improve the VMC calculations is 
through the FN approach)^ where the ground state of 
an auxiliary FN Hamiltonian is obtained. In this case, 
the main approximation relies on the fact that the nodal 
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FIG. 1; (Color online). Energy per hole as a function of the 
doping for J jt = 0.4. Variational (left) and fixed-node (right) 
results are reported for p = and 1 (p = 0, 1 and 2) Lanc- 
zos steps for the wave function with (without) antiferromag- 
netism. The best variational DMRG and iPEPS energies^! 
and the fixed-node with p = 2 are connected by dashed lines 
for a better comparison. 



surface is assigned a priori^ by taking a given guiding 
function that is usually the best variational state. Most 
importantly, the resulting energies are still variational, so 
to have a totally controlled approximation of the original 
problem.* 

In this paper, the guiding function is obtained by op- 
timizing the Jastrow and the mean-field state, with the 
method described in Ref. Then, we find the best 
Lanczos parameters for finally, we perform the 

FN calculations with p = 0, 1, and 2. 

Results. Before showing the results on large systems, 
we would like to mention that a very good accuracy 
on small lattices (where Lanczos diagonalizations can be 
performed) is obtained. We compared our results with 
the exact ones on the 26-site lattice for 2 and 4 holes, and 
different values of J/t (see supplementary material^). 
Both the Lanczos and the FN techniques largely improve 
the variational wave function and the best FN calcu- 
lations (with 2 Lanczos steps) reaches an accuracy of 
{E^^ - E)/Ec^ w 0.002 and w 0.003 for 2 and 4 holes, 
respectively (for J/t ~ 0.4). 

Let us now move to larger sizes and first analyze the 
tendency towards phase separation. In Fig. [U we show 
our results of the energy per hole e{5) — [E{5) — E{0)]/5 
for various cluster sizeS)^ e{6) is a powerful detector for 
phase separation: a monotonic behavior of e{6) vs S indi- 
cates a finite compressibility and a stable uniform phase, 
while a minimum, on finite systems, or a flat behavior in 
the thermodynamic limit, indicate an instability^^ Close 
to half filling, the Pfaffian wave function is considerably 
better than the simple superconducting state, clearly in- 
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FIG. 2: (Color online) Variational results for the variance 
extrapolation on a 162-site cluster, for different numbers of 
holes: p = and 1 (p = 0, 1 and 2) Lanczos steps have been 
performed on the wave function with (without) antiferromag- 
netism. The best fixed-node results are also marked by ar- 
rows. 



dicating a coexistence of pairing and antiferromagnetic 
order liiii^ As the doping increases, the antiferromagnetic 
parameter decreases and eventually vanishes for S ss 0.1. 
The general trend is clear: the increased accuracy of the 
calculation favors the homogeneous state, marked by a 
monotonic behavior of the energy per hole vs the doping. 
In particular, one Lanczos step strongly improves the 
quality of the results, the gain in the FN energy being ap- 
proximately 0.05t, independently of i5. Even the second 
Lanczos step is efficient for these large sizes, providing a 
further energy gain of about 0.02t. We also mention that 
the results obtained with the variance extrapolation are 
consistent with the DMRG and iPEPS onesj^ indeed, 
we have that e{S) = -1.61(1) for 0.03 < S < 0.12. Re- 
markably, we have obtained the same extrapolated values 
(within three error-bars) by using the two wave functions 
with or without antiferromagnetic order, see Fig. [51 How- 
ever, the extrapolated values have too large error-bars 
and cannot be used to study the issue of phase separa- 
tion. 

The application of few Lanczos steps on a given wave 
function is not size consistent; nevertheless, an estima- 
tion of the thermodynamic limit can be attempted by 
considering the largest size, where the p = calculations 
do not show significant size effects. Therefore, we have 
considered p — 2 FN calculations for L = 162 (or even 98 
for 6 > 0.17), which compare well with the best energies 
obtained by DMRG and iPEPS. The latter ones provide 
slightly more accurate energies for 5 ~ 0.1. However, 
considering that all these methods are significantly away 
from the estimated exact energy per hole obtained by 
DMRG and variance extrapolations (i.e., e{6) ~ —1.61), 



FIG. 3: (Color online) Upper panels: local density rii when a 
site-dependent chemical potential with Sjj. — 1.6 [see Eq. Q] 
is added to the variational wave function; the cases with Is = 
12 (a) and 8 (b) are reported. Lower panels: local density 
Hi when a site-dependent potential [see Eq. ((5])] is added to 
the t—J Hamiltonian, with h — 12 and V — 0.2 (c) and 
la — 8 and V = 0.4 (d). Variational and fixed-node results 
are reported for a 12 x 12 cluster and S = 1/8. Insets: the 
difference between the largest and the smallest local density 
(at the fixed-node level) as a function of V. 



this difference looks essentially irrelevant. In contrast 
with DMRG and iPEPS that find a minimum in the en- 
ergy per hole,-^- our best FN approximations do not show 
any tendency to phase separation for any doping, and, 
therefore, represent a thermodynamically stable phase 
corresponding to a well defined variational state. 

Let us now consider the more subtle issue of stripes. 
Recently, DMRG and iPEPS calculations suggested that 
the ground state has charge (and spin) modulations, at 
least close to 6 — 1/8^ Up to now, we have considered a 
uniform mean-field state \^mf), clearly biasing the VMC 
results towards a homogeneous state. Despite the fact 
that the FN method can in principle remove this bias and 
give rise to non-uniform results, we have not found any 
evidence in favor of stripes with this variational ansatz. 

In order to gain some evidence that a charge inho- 
mogeneity is not stabilized in the low-doping regime, we 
add a site-dependent chemical potential in the mean-field 
Hamiltonian 



, 47r 

Mfl. = Mo + oficos [ —X, 



(4) 



where Ri = {xi,yi) is the coordinate of the site i and 
Is equal to 8 or 12. By starting from a finite (5/i, the 
VMC optimization leads to a perfectly uniform state 
with = 0; moreover, FN calculation strongly reduces 
the density modulation present in the original variational 
wave function, see Fig. |31 Although a small inhomogene- 
ity remains in the density profile, the FN energy is always 
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FIG. 4: (Color online) Initial variational ansatz with stripe 
order (a) and fixed-node calculation (b) for charge and spin 
distributions in the 2x8 unit cell of a 16 x 16 lattice. The 
size of the circles and arrows is proportional to the electron 
density and spin along z, respectively. Largest symbols in the 
variational calculations: (n^J = 0.92, (SJ.) = ±0.09. 



higher than the one with (5/i = 0. For these calculations, 
we considered 12 x 12, 16 x 16, and 24 x 24 lattices and 
5 = 1/8. Similar results have been obtained also for 
5 = 1/12 on a 12 x 12 lattice (not shown). 

In order to show the effectiveness and the reliability of 
the FN method to to detect charge inhomogeneities, we 
add a modulated potential directly in the t—J Hamilto- 
nian: 



'47r 

Vr, = y cos ( —X, 



(5) 



Then, we consider a uniform mean-field wave function 
and compute the local density for 12 x 12 and 24 x 24 
lattices and 6 = 1/8. The results are also reported in 
Figs, ini Clearly, the VMC results show a completely 
fiat behavior of the density in different sites; by contrast, 
the FN simulations are able to recover a strongly modu- 
lated density. This fact demonstrates that the presence 
of charge order could be detected by using this approach, 
even when a uniform guiding function is used in the FN 
technique. 

Finally, we can also add a spin structure to the charge 
modulation, so to have: 



(6) 



{riR.) ~ (1 — (5) — (5n cos ^— a;, 
- <5s(-l)«' sin (^^x,^ . 



(7) 



The above structure implies a 2 x unit cell and contains 
the so-called 7r-shift, namely anti-parallel spins across the 
hole-rich sites at Xi = and /s/2. In the following, we 
consider suitable variational parameters inside the mean- 
field Hamiltonian that defines the uncorrelated state (i.e., 
local chemical potentials and local magnetic fields) , such 
to reproduce a stripe with = 8 and take (5 = 1/8 
on a 16 X 16 lattice. Then, we optimize all parameters 
(for each site independently) and observe that the initial 
stripe melts and a perfect uniform state is finally recov- 
ered. Moreover, by performing the FN approach starting 
from a variational state with stripe order, we always ob- 
tain that the charge and spin modulations are reduced 
and a much more uniform state is found, see Fig. |31 we 
also notice that the 7r-shift is replaced by a small defect 
in a weak antiferromagnetic background. 

Conclusions. In this work, we have shown that the FN 
approach is particularly reliable, not only to improve the 
energy of a given variational ansatz, but also to determine 
the density profile of the ground state, in a way that is 
rather independent of the original ansatz. Indeed, the ap- 
proximate FN ground state \'^fn) is not a "brute force" 
variational ansatz, but it represents the ground state of 
a physical Hamiltonian that is different from the exact 
one only in the region where the variational wave func- 
tion is close to zero (namely within the so-called nodal 
region). Operators O that are diagonal in configuration 
space \x) (e.g., related to stripes or antiferromagnetic 
order) are weakly affected by this nodal error. Indeed, 
in the expectation value of O, which takes the form of 
*FAf(^)^2;, the nodal region, where ^'_Fjv(a;) — 0, 
provides a very little contribution, thus explaining the 
reliability of the FN approach. 

We have shown that the FN Monte Carlo, when com- 
bined with few Lanczos steps, is competitive with re- 
cent DMRG and iPEPS calculations, as far as the varia- 
tional energy is concerned. The main outcome is that the 
ground state is homogeneous. No evidence of stripes are 
detected around S = 1/8: at low doping, a uniform state 
is stabilized, containing both superconductivity and an- 
tiferromagnetism. Despite our findings, we have to con- 
clude honestly that the low-doping phase diagram of the 
t—J model is not settled yet, since very accurate methods 
provide very different phases with almost comparable en- 
ergies. We believe that future calculations that employ 
the FN approach on top of iPEPS or DMRG may be 
helpful for the final understanding!^ 

We thank P. Corboz and S.R. White for useful discus- 
sions and for providing us with the DMRG and iPEPS 
results. 
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The energy at half filling has been computed for each clus- 
ter independently. We have seen that this choice minimizes 
finite size effects, that become neglegible for L ~ 162 sites. 
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